Simulating binary neutron stars: dynamics and gravitational waves 



OO 
O 

o 

(N 

C 

c3 



o 
cr 



(N 

> 
o 

(N 
l> 
(N 

OO 

o 
o 



X 



Matthew Anderson, 1 Eric W. Hirschmann, 2 Luis Lehner, 1 Steven L. Liebling, 3 
Patrick M. Motl, 1 David Neilsen, 2 Carlos Palenzuela, 1 and Joel E. Tohlinc 1 

1 Department of Physics and Astronomy, Louisiana State University, Baton Rouge, LA 70803-4001 
B Department of Physics and Astronomy, Brigham Young University, Provo, UT 84602 
'^Department of Physics, Long Island University - C.W. Post Campus, Brookville, NY 11548 

(Dated: February 1, 2008) 

We model two mergers of orbiting binary neutron stars, the first forming a black hole and the 
second a differentially rotating neutron star. We extract gravitational waveforms in the wave zone. 
Comparisons to a post-Newtonian analysis allow us to compute the orbital kinematics, including 
trajectories and orbital eccentricities. We verify our code by evolving single stars and extracting ra- 
dial perturbative modes, which compare very well to results from perturbation theory. The Einstein 
equations are solved in a first order reduction of the generalized harmonic formulation, and the fluid 
equations are solved using a modified convex essentially non-oscillatory method. All calculations are 
done in three spatial dimensions without symmetry assumptions. We use the HAD computational 
infrastructure for distributed adaptive mesh refinement. 



I. INTRODUCTION 

It is widely expected that gravitational waves of suffi- 
ciently strong amplitude will be detected by a new gener- 
ation of gravitational wave interferometers. Binary sys- 
tems composed of compact objects, such as black holes 
and/or neutron stars, are among the strongest expected 
sources of these waves. Advanced gravitational wave de- 
tectors should be sensitive enough to detect the merging 
phase of such binaries. A detailed analysis of the ex- 
pected waveforms from these events will provide valuable 
information not only in the analysis of the received sig- 
nals, but also in the design and tuning of future advanced 
gravitational wave detectors [J Q • 

In relation to these efforts, and beyond the intrinsic 
importance of the two-body problem in general relativ- 
ity (GR), it is significant that recent studies of the bi- 
nary black hole problem have made substantial progress 
in providing waveforms for these mergers (see for in- 
stance 0, H i, S, 0, H). Furthermore, these numer- 
ical results for vacuum spacetimes show a remarkable 
agreement with those obtained with approximation tech- 
niques 0, [13]. This provides considerable support for 
the use of waveforms obtained via approximation tech- 
niques, suitably enhanced by further information from 
numerical simulations, since these can be more easily en- 
coded in a template bank This requires knowing 



the waveforms during the pre-merger, merger and post- 
merger stages and matching them appropriately to obtain 
the continuous wavetrain through the most violent and 
strongly radiative stage of the dynamics. 

For non-vacuum spacetimes, differences in the wave- 
forms may arise from the state of matter describing the 
compact stars, the influence of magnetic fields and re- 
lated phenomena. To fully understand these systems and 
their waveforms, detailed simulations will be required to 
map out the possible phenomenology. 

For the particular case of binary neutron stars in full 
GR, several efforts studying the system in three dimen- 



sional settings have been presented in recent years [EaI131. 

M EE M 

Uj\ , Il8j . However, the complexity and com- 
putational cost of these simulations has permitted inves- 
tigators to consider only a portion of the interesting pa- 
rameter space and several of them have been restricted by 
symmetry considerations. Nevertheless, a number of in- 
teresting problems are beginning to be addressed, includ- 
ing the influence of stiff versus soft equations of state [H[ , 
a possible way to determine the innermost stable circu- 
lar orbit [H| , the dynamics of unequal mass binaries [hsj 
and even the possible existence of critical phenomena in 
the merging system [l7j . 

Further exploration of these systems will require re- 
laxing symmetry considerations, such as axisymmetry or 
equatorial symmetry, and expanding the space of initial 
configurations that can be successfully evolved. More- 
over, the inclusion of additional physics such as magnetic 
fields will be important as these effects may play a major 
role in the resulting dynamics. For instance, the magne- 
torotational instability, which redistributes angular mo- 
mentum in the system, can have a strong influence on 
the multipole structure of the central source and hence 
on the gravitational wave output of the system. 

To date, work on black hole-neutron star binaries has 
been limited to a few cases [ill US [HI- As a result, 
our understanding of this type of system is still in its in- 
fancy. Needless to say, we have even more to understand 
about both types of compact binaries when their environ- 
ments, which may include magnetic fields and radiation 
transport, are included. Indeed, both magnetic fields and 
radiation transport are expected to be key ingredients 
in modeling short, hard gamma-ray burst phenomena 
with compact binaries. Understanding such spectacular 
events requires the addition of these ingredients to the 
computational infrastructure. The resulting numerical 
simulations should allow for new astrophysical insights. 

The present work is intended as the first in a scries 
of studies that examine the evolution of compact binary 
systems in full three dimensional general relativity. To 
this end, we have developed a general computational in- 
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frastructure with solvers for the Einstein and relativistic 
MHD equations that incorporates several novel features, 
which we discuss in the following sections. In Section II 
we describe our formulation of the equations for these 
systems. This includes expressing the Einstein equa- 
tions in terms of a desirable symmetric hyperbolic prop- 
erty [12, [HI and coupling them with the equations of 
relativistic magnetohydrodynamics (MHD) [24j, |25|. Sec- 
tion III presents our numerical implementation, such as 
integration techniques, distributed adaptive mesh refine- 
ment (AMR) , and a tapered grid algorithm that ensures 
stability and considerably reduces spurious reflections off 
artificial internal boundaries [26j]. These ingredients let 
us simulate binary evolutions in which the stars begin 
with wider separations than has been done in earlier stud- 
ies. We can extract gravitational radiation in the wave 
zone and place outer boundaries an order of magnitude 
beyond what has been done previously. As a result, con- 
tamination by boundary effects is negligible. Section IV 
presents a fairly stringent code test by considering the dy- 
namics of a single Tolman-Oppenheimer-Volkoff (TOV) 
solution and extracting the radial oscillation modes of the 
star. Section V describes our main application, namely 
a study of a binary neutron star system without any as- 
sumed symmetries. We follow the dynamics of the system 
from an early non-quasicircular stage to the merger and 
subsequent formation of a neutron star or a black hole. 
We present gravitational wave signals as measured by ob- 
servers placed in the wave-zone and calculated via Weyl 
scalars. Section VI concludes and offers some considera- 
tions for future work. 



II. FORMULATION AND EQUATIONS OF 
MOTION 

The binary neutron star systems considered here are 
governed by both the Einstein equations for the geome- 
try, and the relativistic fluid equations for the matter. We 
write both systems as first order hyperbolic equations. 

In this section we present a brief summary of our for- 
mulation and equations for both the geometry and the 
fluid. More details on our approach to the Einstein equa- 
tions [22j and the relativistic fluid equations 0, [H| can 
be found elsewhere. By way of notation, we use letters 
from the beginning of the alphabet (a, &, c) for spacetime 
indices, while letters from the middle of the alphabet (i, 
j, k) range over spatial components. We adopt geometric 
units where c = G = 1. 



A. Einstein equations 

We write the Einstein equations in a first order reduc- 
tion of the generalized harmonic (GH) formalism. Our 
approach is closely related to the one in [27|, and it 
was used previously in binary boson star evolutions [22l | , 
where additional information can be found. 



We define spacelike hypersurfaces at x° = t = const., 
and define the 3-metric hij on the hypersurfaces. A 
vector normal to the hypersurfaces is given by n a — 
— V a i/| | V a i| |, and coordinates defined on neighboring 
hypersurfaces can be related through the lapse, a, and 
shift, /3\ With these definitions, the spacetime metric 
g a t can then be written as 

ds 2 = g ab dx a dx b (1) 
= -a 2 dt 2 + hij (dx l + j3 l dt) (da? + j3 j dt) . (2) 

Indices on spacetime quantities are raised and lowered 
with the 4-metric, g ab , and its inverse, while the 3-metric 
and its inverse are used to raise and lower indices on 
spatial quantities. 

In the generalized harmonic formulation, the evolved 
variables are 

9ab , Qab = —n c d c g a b , Di a b = dig ab , (3) 

namely the spacetime metric and its temporal and spatial 
derivatives, respectively. Coordinates are specified via 
the generalized harmonic condition 

Ux a = H a (t,x i ), (4) 

where the arbitrary source functions H a (t, x l ) determine 
the coordinate freedom. Although our code allows for 
a general coordinate choice, we choose harmonic coordi- 
nates for the work presented here and set H a (t,x z ) = 0. 
The evolution equations in our GH formalism are 

d t g a b = P k D kab - a Q a b, (5) 
d t Qab = h d k Q ab -ah ij diDj ab 

- a d a H b - a d b H a + 2 a T cab H c 

+ 2ag cd (h ij D ica D jdb - Q ca Q db - g ef T ace T bdf ) 

- ^n c n d Q cd Q ab - a h v D iab Q jc n c 

- 8ira(2T ab - g ab T) 

- 2a a [n a Z b + n b Z a - g ab n c Z c ] 

+ <JiP l {D lab -d igab ), (6) 

d t D iab = (3 k d k D iab - a d l Q ab 

+ -n c n d D icd Q ab + a // " //' /),„ I >,.,,, 

- ci a (Diab - dig a b)- (7) 

Here T ab is the stress-energy tensor and T is its trace, 
T = T a a - Z a is a vector related to the constraints de- 
fined below in Eq. (Ti"Tj) . These variables are not evolved, 
rather they measure the constraint violation and are in- 
cluded in the evolution equations for constraint damping 
purposes [HI]. We also define T abc = g ad T d bc , where T a bc 
are the Christoffel symbols obtained from g ab , given by 

r°6c = ~ 9 ad (Dbdc + D cd b — D dbc ) . (8) 

Note, Di ab are evolved variables in our system, and the 
quantities D 0ab are computed from Q ab and D iab via 

D 0ab = -uQab + l3 k D kab . (9) 
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While the Arnowitt-Deser-Misner (ADM) extrinsic cur- 
vature is not part of the GH system, the fluid equations 
below are written in terms of Kij, which can be calcu- 
lated as 



K, 



1 



1 



(-D(ij)o — (3 k D(ij) k ). 



(10) 



This GH formulation includes a number of constraints 
that must be satisfied for consistency, including the 
Hamiltonian and momentum constraints as well as addi- 
tional constraints that arise in the first order reduction. 
In particular, if we define the four-vector 



2Z a 



pa be rja (± i \ 

-r be g - h [t,x ) , 



(ii) 



it can be shown that the energy and momentum con- 
straints are satisfied if Z a = = dtZ a . The free param- 
eters do and (Ti are chosen to control the damping of the 
four vector Z a (the energy and momentum constraints) 
and the first order constraints, respectively [H, H3|- We 
monitor the Z a during the evolution as an indication of 
the magnitude of the numerical error in the solution. 



B. Perfect fluid equations 

We now briefly introduce the perfect fluid equations. 
Additional information can be found in our previous 
work [24], HH as well as in general review articles [29|, H3] • 

The stress-energy tensor for the perfect fluid is 



Tab = h e u a u b + Pg a b, 



(12) 



where u a is the four velocity of the fluid, h e is the en- 
thalpy, and P is the isotropic pressure. The enthalpy can 
be written 



h e = Po + p a e + P, 



(13) 



where p the rest energy density, and e is the specific 
internal energy density. We introduce the quantities 



W 



-n u a , 



v l = — h\u 3 , 
W J 



(14) 



where W is the Lorentz factor between the fluid frame 
and the fiducial ADM observers and v % is the spatial co- 
ordinate velocity of the fluid. The set of fluid variables 
introduced here are known as the primitive variables, 

W = (p ,v\P) T . 

High resolution shock capturing schemes (HRSC) are 
robust numerical methods for compressible fluid dy- 
namics^ These methods, based on Godunov's seminal 
work -3l|, are fundamentally based on writing the fluid 
equations as integral conservation laws. To this end, we 
introduce conservative variables q = (D, Si, r) T , where 



D = W Po , 
S t = h e W 2 v t , 
t = h e W 2 -P-D. 



(15) 
(16) 
(17) 



In an asymptotically flat spacetime these quantities are 
conserved, and are related to the baryon number, mo- 
mentum, and, in the non-relativistic limit, the kinetic 
energy, respectively. Anticipating the form of the evo- 
lution equations, we also introduce the densitized con- 
served variables 



D = VhD, Si=VhSi, 



Vhr, (18) 



where h = det(/iy). The fluid equations can now be 
written in balance law form 



d t q + d k f\q) = s{q), 



(19) 



where f k are flux functions, and s are source terms. The 
fluid equations in this form are specifically 



d t D + d z 



dtSj 



aD I v l - — 

a 



= 0, 



(20) 



- — ) +VhPh* 



a 3 T 



s., 

\ ]k (SiV k + VhPhi 
-d 3 a{f + D), 



+ S a dj(3 a 



Ka&v 1 + VhKP- 



- S a d a a 
a 



(21) 



(22) 



Here 3 T l j k is the Christoffel symbol associated with the 
3-metric hij , and K is the trace of the extrinsic curvature, 
K = K\. 

Finally, we close the system of fluid equations with an 
equation of state (EOS). We choose the ideal gas EOS 



P=(T-l) Po e, 



(23) 



where T is the constant adiabatic exponent. Nuclear mat- 
ter in neutron stars is relatively stiff, and we set T = 2 
in this work. When the fluid flow is adiabatic, this EOS 
reduces to the well known polytropic EOS 



P 



np 



(24) 



where k is a dimensional constant. We use the polytropic 
EOS only for setting initial data. 



III. NUMERICAL APPROACH 

Our numerical approach to solving the combined equa- 
tions of general relativistic hydrodynamics (GRHD) is 
built upon two extensively tested codes: These were 
written to solve the (1) Einstein equations [12,[23j], and 
the (2) relativistic magnetohydrodynamics (MHD) equa- 
tions [24], [2|| . It should be mentioned that while we do 
solve the full GRMHD equations, in our current work 
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the magnetic field is set to zero. Results with non-trivial 
magnetic fields will be presented elsewhere (32j . 

While both sets of evolution equations are hyperbolic, 
the solutions from each set of equations are quite dif- 
ferent. The Einstein equations are linearly degenerate, 
and therefore we expect smooth solutions to evolve from 
smooth initial data. The fluid equations, on the other 
hand, are genuinely nonlinear, and discontinuous weak 
solutions (shocks) generically evolve from smooth ini- 
tial data [23| . We choose numerical methods adapted 
to the features of each set of equations. The fluid equa- 
tions are evolved with a modified convex essentially non- 
oscillatory (CENO) method, while the Einstein equations 
are evolved using fourth-order accurate difference opera- 
tors that satisfy summation by parts (SBP). These very 
different methods are easily combined by discretizing the 
equations in time using the method of lines. 

We base our code on the HAD computational infrastruc- 
ture for distributed AMR. The Einstein and fluid solvers 
are written in separate modules, which can be used indi- 
vidually or combined. The following sections review our 
methods. 



A. Adaptive mesh refinement using HAD 

The neutron star problem has several important phys- 
ical scales, and each must be adequately resolved to cap- 
ture the relevant dynamics. These scales include (1) the 
individual stars, preferably incorporating some of their 
internal dynamics, (2) the orbital length scale, (3) the 
gravitational wave zone, and (4) the location of outer 
boundaries. In this work, the initial orbital scale is 
on the order of several stellar radii, the gravitational 
waves are extracted at 30, 40 and 50 stellar radii, and 
the outer boundaries of the computational domain are 
placed about 100 stellar radii from the orbital pair to 
reduce boundary contamination of the orbital dynamics 
and gravitational wave signals. The computational de- 
mands required to resolve these different physical scales 
are best met using adaptive mesh refinement. 

We use the publicly available computational infrastruc- 
ture HAD to provide parallel distributed AMR for our 
codes [13, HI] . HAD can solve both hyperbolic and ellip- 
tic equations, and, unlike several other publicly avail- 
able AMR toolkits [H 13, HI US E3> Jt accommo- 
dates both vertex and cell centered algorithms, had 
has a modular design, allowing one to solve different 
sets of equations with the same computational infrastruc- 
ture. Furthermore, solvers for different equations can be 
coupled together, as we have done here with separate 
solvers for the GR and MHD equations, had provides 
Berger-Oliger [4l| style AMR with subcycling in both 
space and time. The HAD clustering algorithm is Berger- 
Rigoutsos [iH, and the load balancing algorithm is the 
least loaded scheme 43] . Refinement in HAD can be trig- 
gered by user-specified criteria, e.g., refining on solution 
features such as gradients or extrema, or refining on trun- 




FIG. 1: The AMR mesh structure at times 0, 84, and 500 
for the pre-merge stage of the simulation with resolution of 
32 points across each star. The simulation had seven levels 
of refinement, five of which are visible here. Simulations were 
performed on 128 processors. 



cation error estimation using a shadow hierarchy. The 
runs presented here use the shadow hierarchy for refine- 
ment, and all dynamic fields are used to estimate the 
truncation error. Some additional fixed refinement re- 
gions are used for gravitational wave extraction in the 
wave zone. As an example, Figure [1] illustrates the re- 
sulting mesh structure at a pre-merge stage in our simu- 
lations. 

had supports arbitrary orders of accuracy [26| , and the 
overall accuracy for the implementation employed here 
is third order for smooth solutions, had implements 
the tapered-grid boundary method for internal bound- 
aries [2fJ. This method is advantageous for two reasons. 
It guarantees stability of the AMR algorithm if the uni- 
grid counterpart is stable as well as significantly reducing 
spurious reflections at interface boundaries. 

Finally, when a fine grid is created during an evolution, 
the geometric variables are interpolated onto the fine grid 
using Lagrangian interpolation. The fluid variables are 
interpolated using weighted essentially non-oscillatory 
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(WENO) interpolation [44]. This interpolation scheme 
is designed for discontinuous functions, and reduces to 
Lagrangian interpolation for smooth functions. 

B. Method of lines 

The numerical methods for the Einstein equations 
(SBP) and the fluid equations (CENO) both specify the 
discretization of the spatial difference operators, giving 
the semi-discrete equations 



Here u represents the set of all variables evolved in both 
the Einstein and fluid equations, and L represents a dis- 
crete spatial difference operator. These ordinary differ- 
ential equations are now discretized in time using the 
method of lines. We choose a third order Runge-Kutta 
scheme that preserves the TVD (Total Variation Dimin- 
ishing) condition [45| to integrate the semi-discrete equa- 
tions 

u« = u n + AtL{u n ), 

M (2) = ^«™ + I m (i) + Ia<L(m( 1 )), (26) 
4 4 4 

3 3 3 v ; 

Using the method of lines for the temporal discretiza- 
tion gives us considerable freedom in choosing numerical 
methods for the spatial derivatives, as well as the ability 
to choose methods of arbitrary orders of accuracy. This 
freedom allows us to naturally and consistently combine 
both the CENO and SBP methods in the GRHD code. 



C. Einstein equations 

As described in our implementation of the Einstein 
equations takes advantage of several techniques tailored 
to the symmetric hyperbolic properties of the generalized 
harmonic formulation we use. At the linear level, these 
techniques guarantee that the full AMR implementation 
is stable. We use second and fourth order spatial deriva- 
tive operators which satisfy summation by parts. These 
operators allow one to obtain a semi-discrete energy esti- 
mate which, together with suitable boundary conditions 
and time integration, ensure the stability of the imple- 
mentation of linear systems (see [Igj], also [13] and ref- 
erences cited therein). Relatedly, we employ a Kreiss- 
Oliger dissipation operator which is consistent with the 
summation by parts property. 

For the outer boundaries, we implement Sommerfeld 
boundary conditions and follow the prescription given 
in 48] . We have also used maximally dissipative bound- 
ary conditions, but found that they led to larger reflec- 
tions at the boundaries which, in turn, corrupt the wave- 
form extraction at late times. 



We set the constraint damping parameters to oo = 
a i = 1. These values were previously used in both binary 
black hole and boson star evolutions, and work similarly 
in the binary neutron star evolutions presented here. For 
the cases discussed here, constraint violations remain un- 
der control during the evolutions. 

Finally, while our GH formalism allows for general co- 
ordinate choices through the source functions H a (t,x' 1 ), 
we adopt H a (t,x l ) — in all the simulations described 
here. Thus, the coordinates adopted are strict harmonic 
coordinates. 



D. Perfect fluid equations 

The perfect fluid equations are integrated using an 
HRSC solver based on the CENO method HH, incor- 
porating some modifications by Del Zanna and Buc- 
ciantini [5Cj . Detailed discussions of our method have 
been presented previously 0, . 

We choose the CENO method to solve the fluid equa- 
tions primarily for two reasons. This means that the dis- 
crete fluid solution corresponds to point values of the so- 
lution and not cell averages. First, it is a finite difference 
or vertex centered scheme. As the Einstein equations are 
discretized with finite differences, coupling these equa- 
tions to the fluid equations with AMR is simplified if both 
sets of variables are defined at the same grid locations. 
Secondly, CENO uses a component-wise decomposition 
(central schemes) of the equations rather than a spec- 
tral decomposition (upwind schemes). Central schemes, 
are more efficient than spectral decomposition schemes. 
Although they are more diffusive at discontinuities, their 
solutions often differ only slightly from those obtained us- 
ing upwind methods. With AMR we can sharply resolve 
all interesting features of the solution. Outflow boundary 
conditions are applied at the physical outer boundary. 

The HLLE flux is used for the numerical flux [5l[ . This 
is a central-upwind method that uses the largest eigenval- 
ues of the Jacobian matrix in each direction. To calculate 
the numerical fluxes, we choose to use piecewise parabolic 
method (PPM) reconstruction for the fluid variables [12] , 
and reconstruct the primitive variables. No dissipation 
or discontinuity detection is used in the reconstruction. 
This is a bit of a departure from the CENO scheme. In 
general, ENO methods use a hierarchical reconstruction, 
where, for example, a second-order reconstruction de- 
pends on an underlying first order reconstruction. We 
have found, at least for the resolutions considered here, 
that CENO often favors a first order reconstruction at the 
center of stars, because of the manner in which candidate 
second order stencils are compared for their similarity to 
the first order reconstruction. This loss of accuracy at the 
center of the star damps the physical quasi-normal oscil- 
lations of the star, and can lead to a long-term growth of 
the central density. PPM gives a superior reconstruction 
for stellar interiors, and therefore we adopt this recon- 
struction here. When the fluid flow is highly relativis- 
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tic, the reconstruction procedure can produce unphysical 
states. When this occurs, we attempt reconstruction us- 
ing a lower order. For example, if PPM fails, then a 
linear minmod reconstruction is attempted, and if this 
fails, then no reconstruction is used. 

A consequence of using HRSC methods is the need to 
go back and forth between primitive, w, and conserva- 
tive, q, variables. While the relation of the conserva- 
tive variables in terms of the primitive variables is alge- 
braic, the transformation that gives the primitive vari- 
ables in terms of the conservative variables is transcen- 
dental. We use a Newton-Raphson solver designed for 
the MHD equations to find the primitive variables [2^ |. 
At grid points where this solver may fail, the primitive 
variables are obtained from neighboring points by linear 
interpolation. The conservative variables are then recal- 
culated at these points from the interpolated primitive 
variables. 

Unphysical states can arise during the evolution of the 
fluid equations. This often occurs in evacuated regions 
of the grid, where truncation errors or effects from finite 
precision arithmetic are significant compared to the fluid 
densities. To compensate for some of these errors, a floor 
is applied to the energy variables D and f as 



D 

f 



max(j), floor), 
max(f, floor). 



(27) 
(28) 



The floor in these runs is set between 1 x 1CP 8 and 5 x 
10~ 9 , which is seven orders of magnitude smaller than 
the central rest mass densities (p c ) of the individual stars. 
The floor value must be small compared to the densities 
in the problem so that the floor does not significantly 
affect the dynamics of interest. Often the effect of the 
floor can only be ascertained by varying it in a series of 
runs. For example, we found that floor values of 10 -7 are 
too large, producing a noticeable increase in p c during the 
evolutions, and changes in the stellar trajectories and the 
emitted waveforms. These errors essentially disappear 
when the floor is 10 -8 , and reducing it further to 5 x 10 -9 
does not change the solutions. Thus, we adopt here a 
floor of 1 x 10~ 8 . 



IV. OSCILLATING MODES OF SINGLE TOV 
STARS 



Mode 


3D GRHD code 
(kHz) 


Perturbation results 
(kHz) 


Relative Difference 

(%) 


F 


14.01 


14.42 


2.88 


HI 


39.59 


39.55 


0.1 


H2 


59.89 


59.16 


1.2 


H3 


76.94 


77.76 


1.1 



TABLE I: Comparison of small radial pulsation frequencies 
for an evolved star using the 3D GRHD code to the linear 
perturbation modes [54|. The polytrope is constructed for 
r = 2 and k = 1. The perturbation results have been appro- 
priately rescaled for k, = 1 55] . The Fourier transform of the 
central density time series is plotted in Figure [3] 



mass of M = 0.14, a circumferential radius R = 0.958, 
and central rest mass density p c = 1.28 x 10 -1 . We evolve 
the data in a dynamic spacetime at different resolutions 
and using different reconstruction methods for the fluid 
variables. Figure [2] shows p c plotted as a function of 
time for three resolutions of 32, 64, 128 points across the 
star. As expected, the oscillations and overall drift in 
p c converge with resolution. This is important both as 
a code test and an indication of the resolution necessary 
to capture some dynamics of stellar interiors. The data 
in Figure [2] were generated using PPM reconstruction. 
We found that first- and second-order CENO reconstruc- 
tions were more diffusive, resulting in larger drifts in p c . 
Consequently, we had difficulty in reproducing the radial 
pulsation modes of the star using these reconstructions. 

To confirm that the code reproduces the expected 
physical behavior, we examine the radial pulsations of 
the star. The modes are calculated from the oscillations 
in p c , and the extracted frequencies are shown in Ta- 
ble HI (Though we present data for the central density 
only here, we have verified that these are global modes 
by examining the time variation of density and veloc- 
ity in the star.) These oscillation modes can be com- 
pared to the known radial perturbation modes [54| . and 
these frequencies are in excellent agreement. Note, to 
make these comparisons we rescale the perturbation re- 
sults as described in the Appendix, which were calculated 
for k = 100. These validations are a stringent test of our 
computational methods and give us considerable confi- 
dence that our code accurately reproduces the physics of 
these systems. 



As a first test of our combined GRHD code we consider 
a single TOV star. Our goal is not only to represent the 
analytic TOV solution, but to accurately reproduce the 
known radial oscillation modes of the star. While the 
TOV solution is spherically symmetric and static, dis- 
cretization effects act as small perturbations that excite 
the normal modes of the star. 

The initial data for this test consist of a T — 2 poly- 
trope with k — 1. (The solution is calculated using a 
modified version of the RNS code of Stergioulas [53j|.) 
The star, in the geometrized units with k = 1, has a 



V. BINARY NEUTRON STARS 

In this paper, we consider two different binary neutron 
star mergers, one resulting in a prompt collapse to a black 
hole and one that results in a differentially rotating neu- 
tron star which persists for a long time (as compared to 
an orbital time close to merger). We evolve the systems 
through several orbits and extract gravitational radiation 
from the orbiting phase and the merger. In the course of 
performing these evolutions, we carefully examine some 
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FIG. 2: This figure shows oscillations in the central rest en- 
ergy density, p c , for a dynamic spacetime evolution of a sin- 
gle TOV star at three different resolutions: 32, 64, and 128 
points across the star. The initial data are for a star of mass 
M = 0.14, circumferential radius 0.958, central rest mass den- 
sity p c = 1.28 x 10 _1 , r = 2, and k = 1. The outer boundary 
of the simulations is 12 stellar radii away from the center of 
the star, and PPM is used to reconstruct the fluid variables. 
While p c increases noticeably for the coarsest resolution run, 
it eventually stabilizes at a higher value, giving a stable con- 
figuration. Results from the Fourier transform of this data are 
given in Table [I] Owing to the computational costs of these 
simulations, the higher resolution runs were not evolved to 
the same end time. In particular, the highest resolution run 
was evolved only until t ~ 65. 




FIG. 3: This figure shows the Fourier transform of the oscilla- 
tions of p c seen in the highest resolution simulation of Figure 
[2] Five distinct peaks are observed; the first four peaks are 
compared with results found via linear perturbation (See Ta- 
ble . The scale of the vertical axis is arbitrary. 
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FIG. 4: This figure shows the maximum value of p in binary 
simulations at four resolutions: 8, 16, 24 and 32 points across 
each star. With fewer than 16 points across the star, the stars 
disperse. For increasing resolutions, the solutions converge. 



numerical questions to ensure the accuracy of our results. 

Initial data for both binaries are set by superposing; the 
initial data for single, boosted TOV stars [5y|. Provided 
that the initial separation between stars is sufficiently 
large, violations in the momentum and Hamiltonian con- 
straints are at or below the truncation error threshold. 
We monitor that this is indeed the case for our chosen 
separations by evaluating the constraints and checking 
that any violations are of the same order as those ob- 
tained for the single stars considered in the previous sec- 
tion. Thus, these data are numerically consistent. The 
boost velocities are smaller than the corresponding Ke- 
plerian velocities for Newtonian circular orbits. Thus, 
our data arc not quasi-circular (as used in and 
they do not correspond to a system resulting from a long, 
slow inspiral. However, these data allow us to both test 
our implementation, as well as to examine how radiative 
effects circularize the orbits. Forthcoming work will con- 
sider initial data taken from post-Newtonian and quasi- 
cquilibrium approaches. 

We extract the gravitational wave information by com- 
puting the Weyl scalar ^4, and for convenience we fur- 
ther decompose r^4 as an expansion in terms of (spin- 
weighted) spherical harmonics 



-2 V 
l,m -Mm- 



(29) 



l.m 



This extraction is done at three different locations from 
the center of mass, and we shift the obtained quanti- 
ties in time to account for the travel time between the 
observers along null rays. These observers are placed 
within the wave zone, and the shift in time is given sim- 
ply by the distance in Minkowski spacetime between the 
observers. As we consider here only equal mass bina- 
ries, corrections for gauge effects should be small [57| . 
An analysis of these effects for different binaries will be 
presented elsewhere |58fl. 

As discussed in [15j . boundary and resolution effects 
can strongly influence the dynamics of these systems. To 




200 400 600 800 

Time 

FIG. 5: The coordinate separation between the stars in a 
merging binary is shown here as a function of time for three 
resolutions. Notice that the merger time, about t — 800, is 
almost the same for the two finer resolutions. 



explore the effects of outer boundaries on the simula- 
tion results, we perform two otherwise identical evolu- 
tions with outer boundaries at different locations. While 
a more detailed discussion of these tests follows below in 
Section IV Al we find that outer boundaries at 80 stel- 
lar radii have negligible influence on the solution. To 
examine resolution effects, we adopt a threshold error 
tolerance for the shadow hierarchy such that the result- 
ing mesh covers each star with a minimum of 16 points. 
While in the previous section we used much higher res- 
olutions to capture the interior dynamics of single stars, 
binary evolutions at similar resolutions here are pro- 
hibitively expensive. Figure [4] gives an indication of the 
minimum resolution required to evolve the binary with- 
out resolving the internal dynamics of individual stars. 
Figure \E\ shows the (coordinate) radial distance between 
the center of the stars versus time for the three different 
resolutions. The trajectories converge as the resolution 
is increased. 



A. Black hole final state 

The first set of initial data gives a binary neutron star 
merger that results in a prompt collapse to a black hole. 
As mentioned previously, the initial data are constructed 
from superposing two equal mass neutron stars with zero 
spin angular momentum. In particular, each star has a 
mass of M = 0.89 Mq, a radius of R = 16.26 km, and a 
central density of 3.24 x 10 14 g/cm 3 . The stars are placed 
initially at the coordinate locations (x, y, z) — (0, ±3, 0) 
with the boost v i = (=F0.08, 0, 0). 

We first investigate possible effects from the outer 
boundaries on the simulation results by performing two 
otherwise identical evolutions with the outer boundaries 
at different locations. In one, the outer boundary is at 
80i?, and in the other it is at 124R. These simulations 
use the shadow hierarchy, and the AMR grid-structure is 
determined by the threshold error parameter. An addi- 
tional set of fixed fine grids is placed at larger distances to 
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FIG. 6: Top Panel. The £ = 2, m = 2 mode of r* 4 extracted 
at 50 stellar radii for binary simulations with different domain 
sizes. The smaller domain is of size ±80 stellar radii while the 
larger is ±124 stellar radii. The two results differ only by a 
small phase and amplitude error which appears late in the 
evolution. For both simulations, the floor value is 1 x 10 -8 . 
Bottom Panel: This includes three plots of the £ = 2, m — 2 
mode of r ^4 extracted at 30, 40, and 50 stellar radii. The 
initial data are described in Section [V Al The domain of the 
calculation is 248 stellar radii across. The signals from differ- 
ent extraction surfaces are shifted in time by the appropriate 
(flat-space) differences between the extraction radii. 



ensure sufficient resolution for computing waveforms. As 
a consequence, the grid-structure in the central region is 
determined dynamically while at far distances it is kept 
fixed. We compare the C*2,2 component of the gravita- 
tional wave signal measured by an observer at a fixed 
coordinate distance, 50i?, for the two computational do- 
mains. These waveforms are shown in Figure [6l which 
shows only small differences in the waveforms at late 
times. Additional tests indicate that these differences 
arise from the location of the exterior, fixed refinement 
boxes. This observation is indicated by the coincidence of 
results obtained with outer boundaries at 100i? and 80i? 
with exactly the same coordinate locations of the exte- 
rior grids. The overall excellent agreement between the 
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FIG. 7: The coordinate trajectory of the center of one of the 
neutron stars as it spirals into a black hole end state. The 
points (filled circles) that have been included along the tra- 
jectory are the coordinate locations of the maximum density. 
These points are shown at intervals of At = 20 in order to 
give an idea of the star's speed. 



wave signals suggests that the influence of the boundary 
location is negligible. 

The dynamics of the subsequent evolution shows a 
clear eccentricity which is reflected both in the gravi- 
tational waveforms (bottom panel of Figure [6]) and the 
coordinate trajectories (Figure [7]). It is worth noting 
that, following the suggestion of [§], a waveform simi- 
lar to Figure [5] can be obtained by using the Newtonian 
quadrupole approximation with the coordinate trajecto- 
ries from Figure [7J In addition, these trajectories are 
similar to those obtained by integrating the 2.5 post- 
Newtonian equations. Finally, as with the black hole 
case reported in the orbital coordinate frequency uj c 
(computed from the coordinate trajectories) is in good 
agreement with the orbital waveform frequency lud (com- 
puted from the dominant mode t — 2, m — 2 of r^), as 
shown in Figure [U 

The eccentricity can be computed using the Newtonian 
definition given in 59] 
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FIG. 8: Orbital frequency of the binary as calculated from 
the numerical evolution in two different ways. u> c is obtained 
by following the coordinate position of the centers of the stars 
while ujd is obtained from the dominant mode of r^4. 



0.2- 



0.1- 



(30) 



where oj p is the orbital frequency at a local maximum 
and uj a the subsequent local minimum. The eccentricity 
of this simulation is shown in Figure [9] To compute this, 
we take each half-cycle and evaluate expression ([3TJ]) . thus 
obtaining a discrete set of values. The first point is clearly 
affected by the initial data adopted, but the subsequent 
points show an overall decrease towards zero. This is ex- 
pected as the gravitational radiation carries away angular 
momentum, and its loss circularizes the orbit. 

Upon merger, the object's pressure and rotation can 
not support the star and it quickly collapses to a black 
hole. As described earlier, our simulations are carried 
out with harmonic slicing which is not singularity avoid- 
ing [60(. Although the lapse collapses to zero as illus- 
trated in Figs. [HI and [TT1 it docs not collapse sufficiently 
fast to avoid numerical problems. As a result, the size of 
the merged object decreases rapidly and the code crashes 
when it can no longer resolve the physical length-scales 
within the allowed maximum refinement levels as shown 
in the final frame of Figure 1101 Ongoing work excises 
a region within a trapped surface to avoid this problem. 
We defer to future work a full analysis of the post-merger 
case and the transition to a quasinormal ringing pattern 
in the radiation 1611. 



B. Differentially rotating neutron star 



200 400 600 800 

Time 

FIG. 9: The eccentricity obtained from Eq. ([30]) . After a 
transient behavior due to the initial configuration, an overall 
monotonically decreasing behavior is seen in the eccentricity 
as the binary orbit becomes tighter. 



In the case where the individual stars are initially sep- 
arated (in coordinate space) by 4i? and boosted with a 
speed of 0.0825 the merger does not give rise to a prompt 
collapse to a black hole, rather it produces a single dif- 
ferentially rotating star (See Figs. [121 4131) . As in the pre- 
vious case, the initial orbital dynamics correspond to an 
eccentric inspiral trajectory. But upon merger, the ob- 
ject's pressure and rotation are sufficient to support a 
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FIG. 10: Snapshots at select times viewed down the z axis 
of the orbiting stars and their subsequent collapse to a black 
hole. These snapshots zoom in on the central region of the 
grid and show only a twentieth of the z = slice of the com- 
putational domain. The stars orbit counterclockwise seven 
times before merging and collapsing to a black hole. The ar- 
rows indicate the fluid velocity. The reference vector in the 
upper right hand corner of each panel has a magnitude of 0.5. 
The color scheme indicates the rest mass density. The plots 
show the simulation at times 620, 760, and 820 as indicated 
in the upper left corner of each image. See Figure [14] for a 
plot of the lapse at the origin as a function of time for this 
system. 



newly formed star. The merged object has a bar-like 
structure that is spinning with a characteristic pattern 
frequency. The real part of the coefficient 6*2,2 of 7^4 
for this evolution (shown here in Figure I15p carries a 
signature of the merger (t approximately from 100 to 
200) and of the resulting spinning bar (t greater than 
250). Qualitatively, the outcome of this evolution agrees 
with the results presented for the fully relativistic sim- 
ulation labeled "E-l" in uM, and even with the results 





FIG. 11: Snapshots of the lapse on the z — plane at times 
400, 600, and 820 for the system presented in Figure [10] The 
contours shown correspond to a = 0.8,0.7,0.6,0.5,0.4, from 
the outermost to the innermost one. At times prior to merger, 
only the first contour value exists. After merger, the lapse 
collapses, indicating the formation of a black hole. Notice the 
essentially circular shape of all the contours except for the 
innermost one at the latest time. The lapse at the origin as 
a function of time is shown in Figure 1141 
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FIG. 12: Snapshots at select times viewed down the z axis of 
the orbiting stars and their subsequent merger into a differen- 
tially rotating star. These snapshots zoom in on the central 
region of the grid and show only a twentieth of the z = slice 
of the computational domain. The stars orbit counterclock- 
wise a couple of times before merging. The arrows indicate 
the fluid velocity. The reference vector in the upper right 
hand corner of each panel has a magnitude of 0.5. The color 
scheme indicates the rest mass density. The plots show the 
simulation at times 43, 116, and 210 as indicated in the upper 
left corner of each image. See Figure[l4]for a plot of the lapse 
at the origin as a function of time for this system. 



from the post-Newtonian SPH simulation labeled "Fl" 
in [62| ; co mpare, for example, our Figure [T51 with Figure 
11 in 12j and Figure 3 in [621 ] . These two earlier simu- 
lations also followed the merger of equal-mass, initially 
irrotational neutron stars having a T = 2 equation of 
state. However, the bar-like structure survives noticeably 
longer in our simulation than in the evolution presented 
in [62j , and in our simulation the radiation signature ap- 
pears to carry more detail about the post-merger dynam- 
ics than in either of these earlier evolutions. Specifically, 
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FIG. 13: Snapshots of the lapse on the z — plane at times 
0, 100, and 200 for the system presented in Figure [12] The 
contours shown correspond to a = 0.9, 0.85 from the outer- 
most to the innermost one. At early times, the contour for the 
lowest value is not present. After merger, though the lapse 
evolves to a slightly lower value, it remains bounded above 
~ 0.75. Notice the essentially circular shape of all the con- 
tours at the latest time. The lapse at the origin as a function 
of time is shown in Figure 1141 
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FIG. 14: The lapse at the origin as a function of time for 
the orbiting polytropes and their merger to either a black 
hole or neutron star. For the sake of comparison, we have 



contact. See Figures [111131 for contour plots of the lapse in 
the two cases. 



the structure discernible in Figure IT51 between the times 
180 and 240 reflects the fact that, as viewed from the co- 
rotating frame of the bar, the bar itself is experiencing 
nontrivial oscillations. 

The neutron star that forms from this merger is 
strongly differentially rotating. In an effort to quantify 
this, in the latter stages of the evolution we fit the inter- 
nal motions of the star to a rotation law of the form, 



fi(r) = 



a. 



1 + Ar 2 sin(6») 5 



(31) 



which has proven to be useful in numerous other in- 
vestigations (see for instance, [H, [H, EH)]). Figure [TBI 
shows the time-dependent behavior of the fitted param- 
eters Vt c and A. We note in particular that the ratio of 
the central and near-surface value of f2 at the equator is 



0.34. 



FIG. 15: The merger waveform for the collision resulting in a 
single compact star extracted at three different stellar radii: 
30, 40, and 50. The domain of the simulation is ±152 stellar 
radii. After the merger a transient behavior is observed. In 
particular, the features at t ~ 180, 240 result from marked os- 
cillations in the produced bar-like configuration (as seen in the 
co-rotating frame). Afterwards the gravitational waves due to 
the spinning bar exhibit a clear frequency at ~ 12.8 kHz. 
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FIG. 16: The fitted values of £l c and A as determined from the 
fluid's tangential velocity. The merger takes place at about 
t ~ 140 after which the angular velocity rises during a tran- 
sient stage and then slowly decreases. 



VI. CONCLUSION 



Neutron stars will be important sources of gravita- 
tional waves for the next generation of gravitational wave 
detectors. While waveforms from neutron star binaries 
are weaker than those produced by binary black holes 
due to the allowed neutron star masses, their signals are 
expected to be richer, as the gravitational waves will also 
carry information about the matter. Indeed, gravita- 
tional waves are expected to become an important probe 
of neutron star physics, addressing questions such as the 
equation of state for nuclear matter and the nature of 
progenitors for short, hard gamma-ray bursts. 

We have constructed a code that solves the combined 
Einstein and fluid equations in three spatial dimensions, 



with no symmetry assumptions, and we use HAD for 
distributed AMR. AMR is an essential element of our 
method, as it allows us to place the outer boundaries far 
from the binary, while the shadow hierarchy allows us 
to refine each star individually without a priori assump- 
tions about their motion. We have carefully verified our 
numerical results by performing runs at different reso- 
lutions, using grids with different physical outer bound- 
aries, extracting ^4 at different radii, and varying the 
floor applied to the fluid densities. Moreover, we studied 
the radial pulsation frequencies for a T = 2 polytropic 
TOV star, finding excellent agreement between our re- 
sults and the expected perturbative values. The success- 
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ful conclusion of these tests gives us confidence in the 
physical results obtained from our code. 

As a first application of this code in a demanding sce- 
nario, we present a detailed study of two binary neutron 
star mergers, one resulting in a final black hole and the 
other a final neutron star. In both cases we examine 
the gravitational wave emission by extracting the I = 2, 
m = 2 mode of r^. ^4 is extracted sufficiently far from 
the binary within the wave zone, and extraction is done 
at three different radii. In the first case, '5 4 is extracted 
up until the lapse collapses, and in the second case the 
wave signal is extracted until a final differentially rotat- 
ing star is reached. A comparison to a post-Newtonian 
analysis allows us to understand better the gravitational 
wave signals and the orbital kinematics, such as orbital 
trajectories, frequencies, and eccentricities. For example, 
the initial data describes an eccentric orbit. The effect of 
the eccentricity can be observed in the alternating pat- 
tern of larger and smaller extrema in '5 4 as well as a 
modulation in the observed wavelength. Both features 
are expected from a Post-Newtonian analysis of an ec- 
centric orbit. The orbits circularize through gravitational 
wave emission, and the solution around the time of col- 
lapse is largely spherically symmetric. In the second case, 
the neutron star merger results in a large strongly differ- 
entially rotating star. The observed maximum density 
after the merger does not lie at the origin but oscillates, 
in the co-rotating frame, in a bar-like fashion in between 
~ 0.2i?fi na i and ~ 0.4Rfi na i (with i?fi n ai the equatorial 
radius of the merged object). 

The work presented here raises additional questions 
that we will pursue in a continuing research program. 
For example, we will continue to study the ringdown of 
the final black hole formed in the first merger. Studies of 
the differentially rotating star formed in the second case 
are continuing to determine whether this star eventually 
collapses to form a black hole. We will also examine a 
broader class of initial data, including quasicircular and 
unequal mass binaries. As mentioned previously, we also 
are investigating the effect of magnetic fields on the mas- 
sive compact object formed in a merger and its possible 
subsequent collapse. These results will be published in 
subsequent papers. 



Appendix 

It is customary in general relativity to adopt ge- 
ometrized units G = c = 1, such that all quantities, in- 
cluding mass (M) and time (T), have units of length (L). 
Vacuum solutions are invariant under changes in this 
fundamental length scale L. A quantity X that scales 
as L ; M m T* can be converted into geometrized units by 
multiplying with the factor c* (G/c 2 ) m . After the con- 
version to geometrized units, X scales as L l+m+t . Most 
equations of state break this intrinsic scale-invariance, 
and the fundamental length-scale must be fixed by ad- 
ditional choices. Once the new scale is chosen, trans- 



formations between geometrized and physical units can 
be easily made. In the following, we summarize the ba- 
sic procedure detailed in [55[ to account for the proper 
scaling of quantities. 

The polytropic EOS ([24]) is specified by the constants 
{k, F}, and the quantities obtained when using a partic- 
ular set {kijTi} can be scaled to those obtained using a 
second set {^2^2} by the factor 



Li _ ^V 2 ^- 1 ) 
L 2 ~ K2 i/2(r 2 -i) 



(A-l) 



There are two common approaches in the literature to set 
this additional length scale. The first one is obtained by 
fixing a constant physical quantity, e.g., the solar mass 
Mq = 1, and from it deduce the appropriate conver- 
sion factors. That is, if a quantity X has dimensions of 
L l M m T t , its dimensionless counterpart, X, is obtained 
from the following equation: 



X 



G Mq 



l + L 



M r r 



-X. 



(A-2) 



There is still the freedom to choose k, and all dimensions 
are scaled with this parameter. Usually the choice k = 
100 is preferred because it leads to physical units which 
are close to the current observations. For instance, TOV 
stars constructed with these parameters have a maximum 
stable mass of M max = 1.64Mq with a radius of -R max 
1 4 . 1 I km. 

The second method for choosing the length scale is ex- 
plained in detail in (55l |. and is more involved. It is based 
on fixing the maximum stable mass for a family of solu- 
tions (with given {k = l,r}) to a physically motivated 
value. Thus, a quantity X with dimensions L z M m T* is 
obtained by using the relation: 



X = k x c y G z X, 



(A-3) 



where 



l + m + t 
2(T - 1) ' 
I + 3m + t 



(T - 2)1 + (3T - 4)m - t 



r - 1 



(A-4) 



In this method k has dimensions. We now identify the 
maximum stable mass for the given polytrope to some 
physical maximum mass. For a neutron star, the ob- 
served maximum mass is M max = IAMq. 

Although this second method for fixing the funda- 
mental length scale generally leads to different results 
from the first, it can be checked that for T = 2 both 
methods (the first one with k = 100, while the second 
one always has K = 1) provide the same scaling fac- 
tors when the physical maximum stable mass is set to 
M = 1.64M0. Since the dimensionless maximum stable 
mass is M = 0.164, Eq. (|A-3|) can be solved for k with 
{I = 0,m = 1, t = 0}, giving k = 1.456 x 10 5 cm 5 / (gs 2 ). 
With this value, ([0)1 can again be used to recover the 
dimensions of any quantity. 
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